Next-to-Leading Order QCD Corrections to Heavy Quark Correlations in 
Longitudinally Polarized Hadron-Hadron Collisions 
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We present a comprehensive phenomenological study of heavy flavor distributions and correlations 
in longitudinally polarized proton-proton collisions at BNL-RHIC. All results are obtained with a 
flexible parton-level Monte Carlo program at next-to-leading order accuracy and include the frag- 
mentation into heavy mesons, their subsequent semi-leptonic decays, and experimental cuts. Next- 
to-leading order QCD corrections are found to be significant for both cross sections and double-spin 
asymmetries. The sensitivity of heavy flavor measurements at BNL-RHIC to the gluon polariza- 
tion of the nucleon is assessed. Electron-muon and muon-muon correlations turn out to be the 
most promising observables. Theoretical uncertainties are estimated by varying renormalization 
and factorization scales, heavy quark masses, and fragmentation parameters. 
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I. MOTIVATION AND INTRODUCTION 

Recent results from longitudinally polarized lepton- 
nucleon scattering experiments [l|, [2, El, 3 and, in par- 
ticular, for single-inclusive pion and jet production in 
helicity-dependent proton-proton (pp) collisions at the 
Relativistic Heavy Ion Collider (RHIC) [5j,l6| have started 
to put significant limits on the amount of gluon polariza- 
tion in the nucleon @, [|| . 

This is best quantified in a "global QCD analysis", 
which treats all available experimental probes simulta- 
neously and consistently at a given order in the strong 
coupling a s in perturbative QCD (pQCD). It allows one 
to extract the set of universal, spin-dependent parton 
distribution functions, defined as 



is given by 



(1) 



that yields the optimum theoretical description of the 
combined data. In (JTJ) , / + (/_) denotes the probability 
of finding a parton of flavor / = q,q, g at a resolution 
scale /i with light-cone momentum fraction x and helicity 
+ (— ) in a proton with helicity +. 

Our current understanding of the spin structure of the 
nucleon is derived from Eq. ([1]) by taking the first mo- 
ments of the densities Af(x,p). These quantities enter 
the helicity sum rule of the nucleon along with the con- 
tributions from the orbital angular momenta of quarks 
and gluons [7|. Specifically, the total gluon polarization 
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Ag(x, p)dx 



(2) 



and the challenge is to precisely map the gluon helicity 
density Ag(x, //) in a wide range of x in order to minimize 
extrapolation uncertainties in the first moment Ag(fi). 

A first global QCD analysis of polarized parton den- 
sities Af(x,p,) at next-to-leading order (NLO) accuracy 
was completed recently [8j. It was based on the world- 
data on polarized inclusive and semi-inclusive deep- 
inelastic scattering, which are pivotal in constraining the 
quark and antiquark densities [3] , as well as on the latest 
RHIC pp measurements [fl @ mentioned above. The con- 
clusion is that available results from all lepton-nucleon 
scattering experiments and the RHIC spin program are 
in nice agreement. This underpins the notion of fac- 
torization also for spin-dependent hard scattering pro- 
cesses, which is the foundation for most pQCD calcula- 
tions and their predictive power. The polarized gluon 
density Ag(x, /j) turns out to be compatible with zero in 
the range of momentum fractions, 0.05 < x < 0.2, acces- 
sible to experiments so far. However, it is still impossible 
to give a reliable estimate for the total gluon polarization 
Ag(fi) [8j. A significant contribution to the integral in 
Eq. ([2]) can still come from the unexplored small x region. 
Hence, the fundamental question of what constitutes the 
proton spin still remains largely unanswered, despite the 
fact that impressive progress, both theoretically and ex- 
perimentally, was made in the past two decades. Back 
then, it was discovered that only an unexpectedly small 
fraction, about a quarter, of the proton's spin can be at- 
tributed to the intrinsic spin of quarks and antiquarks 

0- 

Narrowing down the uncertainties on Ag(x, fi) and, at 
the same time, extending the range in x continues to be 
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the main objective of experimental efforts in the years 
to come, utilizing both longitudinally polarized lepton- 
nucleon and proton-proton scattering. With higher lu- 
minosities becoming available at RHIC, less inclusive 
final-states like jet-jet correlations will be instrumental in 
achieving this goal as they give a much better handle on 
the x range probed in experiment Q. Also, rare probes 
like prompt photons and heavy quarks come into focus. 
Both will not be able to compete with single-inclusive 
pion or jet measurements with respect to statistical pre- 
cision, but they follow rather different underlying QCD 
hard scattering dynamics. Therefore, such measurements 
are crucial for further testing and establishing the univer- 
sality of helicity-dependent parton densities and hence for 
our understanding of the spin structure of the nucleon 
and QCD in general. 

In this paper, we present a comprehensive phenomeno- 
logical analysis of open heavy flavor production in lon- 
gitudinally polarized pp collisions at RHIC. Until ver- 
tex detector upgrades are in place, RHIC experiments 
will identify heavy quarks through their semi-leptonic de- 
cay electron or muon spectra, which receive contributions 
from both charm and bottom hadron decays, or by direct 
reconstructions of hadronic D meson decays. Since heavy 
flavors are a versatile probe of high-density medium ef- 
fects in nucleus-nucleus collisions, such as modifications 
of the transverse momentum spectra [To| , various refer- 
ence data have been taken at RHIC in unpolarized pp 
collisions (lH [l2j]. Similar measurements are intended 
with longitudinally polarized protons [l3j . 

To reduce the uncertainties from deconvoluting ex- 
perimental results for decay lepton spectra back to the 
heavy quark level, all theoretical calculations should be 
done as close as possible to the observational level. This 
was achieved, e.g., in a recent phenomenological study 
of unpolarized charm and bottom production at RHIC 
[T3 |. For the simplest example of a single-inclusive elec- 
tron spectrum from semi-leptonic decays of a heavy me- 
son Hq, the corresponding invariant cross section takes 
schematically the following form 

Ej^=E Q d ^^®D^®f»^, (3) 
dpi Pq 

where the symbol <g> denotes a convolution. The cross 
section d(A)a^ for the production of a heavy quark Q 
with mass toq, energy Eq, and momentum pq in (polar- 
ized) pp collisions can be evaluated within pQCD. NLO 
QCD corrections, which are essential for any meaningful, 
quantitative analysis, are known, both in the unpolar- 
ized fl5l[l6| and polarized [l7| case, for quite some time. 
We note that the longitudinally polarized hadronic cross 
section is defined as the combination 

dAa^ = ^[da^ + -da^_], (4) 

where the subscripts ± label the helicity states of the 
colliding hadrons. The result for d(A)a^ depends on 



the choice of non-perturbative (helicity-dependent) par- 
ton densities, the value of ttiq, and on the unphysical 
factorization (fJ,f) and renormalization (/i r ) scales. The 
sensitivity of the cross section to variations of fj,f, r can be 
taken as a rough estimate of the theoretical uncertainty 
due to the truncation of the perturbative series at a cer- 
tain order. Likewise, variations of uiq contribute to the 
theoretical ambiguities as well. We will assess all these 
sources of uncertainties in our detailed numerical studies. 

The other two ingredients to Eq. ([3]), are the fragmen- 
tation D®^ Hq of the heavy quark Q into a heavy me- 
son Hq and the semi-leptonic decay f H Q~^ e of Hq into 
the experimentally observed electrons e. Since uiq cuts 
off final-state collinear singularities associated with the 
heavy quark, its hadronization D®^ H Q is fundamentally 
different from those for light quarks and gluons. In the 
latter case, scale-dependent parton-to-hadron fragmenta- 
tion functions have to be introduced by virtue of the 
factorization theorem. The non-perturbative transition 
Q — > Hq is described by various phenomenological mod- 
els for a scale independent function D®^ Hq , whose pa- 
rameters are determined from fits to e + e~ data [l9| . For 
our phenomenological studies, we use the functional form 
proposed in Ref. [20| . with its single parameter taken in 
the range given in [19(. In addition, a fixed order pQCD 
calculation of d(A)a^ in Eq. ([3]) can be supplemented 
by all-order resummations of quasi-collinear logarithms 
of the form a™ \og n (p® / itiq) [2lj], which can be large 
if the transverse momentum p^ of the produced heavy 
quark is much larger than its mass. For the time being, 
we do not pursue similar resummations for the polarized 
hadroproduction of heavy quarks since p® ~ tuq for all 
phenomenologically relevant applications at RHIC. We 
adopt the parameterization of the semi-leptonic decay 
spectrum f H Q^ e obtained in Ref. [3] from a fit to BaBar 
and CLEO data [22j and used in unpolarized analyses of 
heavy quark production at RHIC [14| . 

The expression in Eq. ([3]) can be easily generalized to 
the important case of heavy quark or decay lepton cor- 
relations. As will be demonstrated in some detail be- 
low, such measurements appear to be more promising 
for accessing the gluon polarization at RHIC than single- 
inclusive decay electron or muon spectra. To make theo- 
retical calculations for such observables feasible at NLO 
accuracy, we develop a flexible parton-level Monte Carlo 
program to perform all phase-space integrations numer- 
ically. For the subtraction of soft and collinear diver- 
gences present at intermediate stages, we follow closely 
the methods devised in Ref. [23[ for the computation of 
heavy flavor correlations in unpolarized hadron-hadron 
collisions. 

Our Monte Carlo code is capable of computing any in- 
frared safe heavy flavor cross section at O(af), including 
correlations of the QQ pair and control of the accompa- 
nying jet, with the same kinematic cuts as used in exper- 
iment. The hadronization of the QQ pair into heavy 
mesons and their subsequent semi-leptonic decays are 
modeled as outlined above. Our results complement and 
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significantly extend previous spin-dependent NLO calcu- 
lations of single-inclusive heavy (anti)quark yields and 
of the heavy quark charge asymmetry based on largely 



although not pursued in this work, with the associated 
hard jet present for the first time at O(af). This signif- 



analytical methods [171 |24| ]. where any information on 



the partonic recoil system was lost. The required spin- 
dependent matrix elements squared at 0(a^) for produc- 
ing a QQ pair plus a light parton, 

99 -> QQg , qq -> QQg , gq(q) -> QQq{q) , (5) 

as well as the virtual corrections to the 0{a 2 s ) Born cross 
sections, 



99 



QQ , qq-^QQ ■ 



(0) 



can be taken from Ref. [17|. Expressions obtained in 
d = 4 — 2s dimensional regularization are required only 
in the singular regions of phase-space, and e can be set 
to zero otherwise. 

After appropriate modifications, the results obtained 
in this work can be used also as the "resolved" pho- 
ton contribution to the spin-dependent photoproduction 
of heavy quarks at NLO. A similar parton-level Monte 
Carlo program including both direct and resolved photon 
processes will be presented in a forthcoming publication 
25]. It will allow one to include existing data on spin- 
dependent charm photoproduction Q into future global 
QCD analyses of polarized parton densities. NLO ex- 
pressions for the point-like, "direct" photon part of the 
cross section can be adapted from (26[; see also [27j . 

The outline of the paper is as follows: in Sec. II we 
briefly review some of the technical aspects of setting 
up a parton-level Monte Carlo program for heavy flavor 
production in polarized hadron-hadron collisions at NLO 
accuracy. Some additional technical details can be found 
in the Appendix. In Sec. Ill we present a detailed phe- 
nomenological study of heavy flavor production in po- 
larized pp collisions at RHIC, focussing on experimen- 
tally relevant decay lepton, electron or muon, distribu- 
tions and correlations within the kinematic acceptance 
of the PHENIX and STAR experiments at RHIC. We as- 
sess theoretical uncertainties due to variations of /// and 
fi r , the heavy quark mass ttlq, and parameters related 
to modeling the fragmentation process. Finally, we give 
predictions for double-spin asymmetries and discuss their 
potential in further constraining helicity-dependent par- 
ton densities, in particular, Ag(x,fj,). We summarize our 
results in Sec. IV. 



II. TECHNICAL FRAMEWORK 

We consider heavy quark hadroproduction in longi- 
tudinally polarized pp collisions at O(af) in QCD. All 
phase-space integrations are performed numerically with 
Monte Carlo techniques. This enables us to compute any 
observable involving heavy quarks within experimental 
acceptance cuts, including single-inclusive distributions, 
correlations among the heavy quark-antiquark pair, and, 



icantly extends available calculations [17|, [24( based on 
largely analytical integrations over the variables charac- 
terizing the partons recoiling from the observed single- 
inclusive heavy quark or antiquark. 

In the integration of the fully exclusive partonic cross 
sections for the processes in Eqs. and ^ at O(o%), 
one has to deal with ultraviolet, infrared, and collinear di- 
vergences, which have to be eliminated before any numer- 
ical approach can be used. To this end, we follow closely 
the subtraction method devised and used in Ref. [23j to 
compute the unpolarized production of a QQ pair. The 
subtraction method is based on adding and subtracting 
counter terms which approximate the real emission pro- 
cesses in (|S|) in the singular regions of phase-space and 
are integrable with respect to the momentum of an un- 
resolved parton. In this Section, we briefly review the 
technical aspects relevant to extend and apply the sub- 
traction method of [23| to heavy quark hadroproduction 
in longitudinally polarized pp collisions and discuss the 
numerical implementation. For further details, we refer 
the reader to Ref. [23[ ■ We note that a general formula- 
tion of the dipole subtraction method for NLO calcula- 
tions with massive partons in QCD and supersymmetric 
QCD has been developed in [281 ]. 

Assuming, as usual, factorization, the cross section ^ 
for producing a heavy QQ pair in longitudinally polarized 
pp collisions at a center-of-mass system (c.m.s.) energy 
y/S can be written as a convolution, 

dAa Q = ^2 J dxidx 2 Af a (xi,pf)Af b (x 2 ,Pf)S 

a,b 

x dAa a b(xi,x 2 ,S, mQ,ki,k 2 ,Hf,^r) , (7) 

where the Af a ^(x a ,b, M/) denote the spin-dependent par- 
ton distribution functions of flavor a, b at momentum 
fraction x a ,b and scale pf, as defined in Eq. {l}. The 
sum in ([7]) is over all contributing partonic processes 
ab — > QQc to 0(al) with dAa a b the associated polar- 
ized hard scattering cross sections. They are defined in 
complete analogy to Eq. ([4]) and can be computed per- 
turbatively as a series in the strong coupling a s . Parton 
c is either a gluon or a light (anti)quark producing the as- 
sociated jet possible at 0{a\). k\_ 2 denote the momenta 
of the heavy quark Q and antiquark Q with mass rriQ, 
i.e., k 12 = itiq. 

The required spin-dependent matrix elements squared 
at C(a^) in d = 4 — 2e dimensional regularization for 
the processes in ([5]) and ((6|) can be taken from Ref. [re- 
starting from the NLO level, the subprocess cross sec- 
tions dAa a b in ([7|) depend explicitly on the renormal- 
ization and factorization scale p r and fj,f, arising from 
the subtraction of ultraviolet and collinear singularities, 
respectively. Infrared (soft gluon) divergences cancel 
among real emission and virtual loop corrections. 

In Eq. 0, S is the "measurement function" used to 
define the observable one is interested in. One can think 
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of S as being a set of step functions implementing the ex- 
perimental cuts imposed on the final-state particles and 
selecting a certain bin in a histogram. As mentioned 
in the Introduction, charm and bottom quarks are cur- 
rently detected only indirectly at RHIC, mainly through 
the semi-leptonic decays of the produced heavy D and B 
mesons. Thus, the cross section ([7]) at the heavy quark- 
level is not yet sufficient for comparing theory with ex- 
perimental results. As indicated in Eq. ([3]), one needs to 
convolute the parton-level results for dAa® with addi- 
tional phenomenological functions D^"* Hq and f H Q^ e 
describing the hadronization into a heavy meson Hq and 
the semi-leptonic decay of Hq into the observed lepton, 
respectively. Our flexible parton-level Monte Carlo pro- 
gram not only performs the phase-space integrations for 
arbitrary S for any infrared safe observable but can also 
account for the semi-leptonic decays of the heavy quark 
pair into electrons and muons. We specify our choice for 
D Q^H Q and jH Q -,e in Scc. lTrTA"! 

For the implementation of Eq. in a numerically 
efficient integration it is convenient to express the three- 
body phase-space and the matrix elements squared for 
the 2^3 processes listed in in terms of variables 
in which soft and collinear singularities can be identified 
easily. Instead of choosing the usual set of five indepen- 
dent scalar products (or Mandelstam variables) of the 
parton momenta in ab — ► QQc, this is achieved by in- 
troducing x, y, 0i, 8 2 , and s — X1X2S [231 ] . They are 
defined as follows: x = (fei + k^) 2 /s, the invariant mass 
of the QQ pair scaled by the available partonic c.m.s. 
energy squared, i.e., p = Atuq/s < x < 1, and y is the 
cosine of the angle between the z-direction, aligned with 
the spatial direction of parton a, and k%, the momentum 
of parton c, in the c.m.s. of the incoming partons, i.e., 
— 1 < ^/ < 1. Soft and collinear regions of phase-space are 
associated with x = 1 and y = ±1, respectively. Both 
81 and 02 do not matter for this discussion. They range 
between and n and are used to parameterize the spatial 
orientation of fcj j with respect to the plane span by the 
other three momenta in the c.m.s. of the QQ pair, see 
[23| for an explicit parameterization of the momenta. 

The d-dimensional three-body phase-space expressed 
in terms of the variables x, y, 6\, 62, and s reads 



1 



2 -9+6e 7r -4+2e s l-2 £/3 l~2e 



r(i-2 £ ) 

x x-^l-xf-^dxil-y^dy 



x sin 



1— 2s , 



'i d6\ sin 



-2e , 



•d6 2 



(8) 



which agrees with the result in [23j and where we have 
introduced /3 X = [1 — 4to 2 ->/(sx)] 1 / 2 . T(z) represents the 
Gamma function. 

The contribution of the 2^3 real emission processes 
in (0 is then given by 



dAa ab = A|M ab | 2 dPS 3 , 



(9) 



where the spin-dependent amplitude squared, A|M a b| 2 , 
includes the partonic flux factor l/(2s) and is summed 



over final-state color and spin degrees of freedom and 
averaged over the color of the interacting partons a, b 
[13]. Soft (x — 1) and collinear (y = ±1) singularities in 
A|M a fc| 2 appear as 



A|M„ 



Af ab (s,m Q ,x,y, 6>i6> 2 ) 
s 2 (l-x) 2 (l-y 2 ) 



(10) 



where Af ab is regular for x — 1 and y = ±1. The qg pro- 
cess in ([5]) can have only collinear singularities at NLO. 
Due to the finite mass tuq, there can be no collinear 
gluon radiation from a heavy quark ("dead cone"). 

Upon inserting ((8]) and (fT0|) into |9]), one can proceed 
by expanding the resulting (1 — a;) _1_2e and (1 — y 2 )^ 1 ^ 6 
for small e as shown in Ref. [231 ]. 



(1-x) 



-l-2e 



-4£ 



-5(1 - x) + 



2t 



2s 

log(l - x) 
1-x 



1 



1-x 



+ 0(e 2 ) 



(i-y 2 ) 



2\-l- 



.[S(l+y)+S(l-y)] 



(2uy 

2e 



1 



1-2/ 



1 



1 + y 



0(e), 
(11) 



where (3 — y/1 — p. Explicit expressions for the distri- 
butions [1/(1 - x)]p, ... in HI]) are collected in Eq. (J30j) 
of the Appendix. The choice of the parameters p and lu 
is to some extent arbitrary and will be discussed at the 
end of this Section. Using Eq. (fTTjl . the subprocess cross 
sections for ab — > QQc at 0(a? s ) can be decomposed as 



dAa a b = dAa 



(b) 



dAa 



dA, 
dAa 



j.(c+) 

ab 



+ dA 
dAa 



ab 

(/) 



fl2) 



Here, dAa^ and dAa [ a " b ' denote the 0(aj) Born con- 
tribution and the O(a^) one-loop corrections to the gg 
and qq scattering processes in (j6]), respectively. Analytic 
expressions for the virtual contributions in d dimensions, 
with ultraviolet divergences being subtracted at a renor- 
malization scale fi r , have been obtained in Ref. [13 ] . 

In Eq. (fT2|) . dAa^j is the soft component of the gg 
or qq scattering cross section, which can be either eval- 
uated by explicitly taking the soft gluon limit of the 
full d-dimensional matrix elements squared computed in 
Ref. [13] or constructed using general properties of soft 
gluon emission in QCD, see, e.g., [23]. In the limit x — ► 1, 
the kinematics simplifies, and phase-space integrations 
can be performed analytically. The relevant integrals are 
the same as for unpolarized heavy flavor hadroproduc- 
tion and can be found, e.g., in App. A of Ref. [23§. One 
obtains 



dAa 



M 

ab 



T(l 



r(i - 2s) 



-4c 



Af^(s,m Q ,e 1 )dPS2. 



(13) 
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Explicit expressions for dAa^ and, for completeness, the 
standard two-body phase-space factor dPS2 in d dimen- 
sions are given in the Appendix. 

All 2^3 processes in (|5|) exhibit singularities related 
to collinear splittings off the incoming partons. Again, 
for such configurations the kinematics collapses to the 
simpler case of 2 — > 2 scattering, and these contribu- 

i.(c±) 



where 



tions, summarized by dAa ab in Eq. (|12p . can be eval- 
uated analytically. As for the soft contribution, one can 
either start by taking the collinear (y — + ±1) limit of the 
full, rf-dimensional 2 — > 3 matrix elements taken from 
Ref. [ItJ or by deriving the expressions from scratch. Af- 
ter combining the relevant matrix elements with g?PS3 in 
Eq. ([8]), taking the limit y — ► ±1, and integrating over di 
one obtains 



dAa 



-(47r)^r[l+ £ ](- 



4e 



-dPS? 



1 - x 



2s 



log(l 



1 -x 



^fab i s ' m Q' x ' 



(14) 

where dPS% — dPS2\ s _ txg dx. The superscript ± in 

dAa^^ distinguishes the two configurations with y = 
+1 and y = — 1, where parton c is emitted collinear ly 
to the momentum of parton a and 6, respectively. The 
relevant A/^^ are again collected in the Appendix. 

The last term in Eq. ([12)) . dAa^ b \ contains all the 
finite contributions after using the expansions (|lip for 
(1 — a;) _1_2e and (1 — y 2 )~ 1_e , and the phase-space inte- 
gration can be performed numerically in four dimensions, 
i.e., with e — ► 0. One obtains 



dAa 



(/) 

ab 



l 



1-x 

x sin 9\d9i d0 2 dx dy 

x Af ab (s,m Q ,x,y, 9i,0 2 ) ■ (15) 

As can be seen, all soft and collinear singularities are 
regulated by the p- and ^-prescriptions defined in (|3H)) . 

The resulting 1/e divergence in (Ti"4"|) assumes the form 
dictated by the factorization theorem, i.e., a convolution 
of d-dimensional helicity-dependent LO splitting func- 
tions APij(x) and Born matrix elements A|Af Q {,| 2 . Due 
to the collinear splitting, the latter have to be evaluated 
at a "shifted kinematics" where parton a (or b) carries 
only a fraction x of its original momentum, i.e., s — ► xs 
and dPS 2 -> dPS 2 ; see Eqs. (|55 ]l -(|i5 |l in the Appendix, 
where, for convenience, also the Born cross sections and 
the LO APij(x) are listed. Collinear singularities are fac- 
torized into the bare parton distribution functions at a 
scale \if by adding an appropriate "counter cross section" 
to (|14p which to O(al) schematically reads 



2tt 



AVia (X, fMf)dA(T$ (xs 



+ AV lb (x,fi f )dAa { ^(xs) 



V tj (x,fi f ) = AP l3 (x) [- - + lE - In 4tt + In -£] + Ag zj (x) . 

(17) 



(16) 



The sum in (|16|) is over all possible collinear configura- 
tions involving one of the initial-state partons a, b. The 
argument xs of the Born cross sections dAa^ in (|16|) 
indicates that they have to be evaluated at the shifted 
kinematics as discussed above. In Eqs. (fTB"]) and (|17p . 
the Euler constant 7^ and ln47r, both, like the scale 
H, artifacts of dimensional regularization, are subtracted 
along with the 1 / e singularity. The factorization scheme 
is fully determined by the choice of Agij, for which we 
take Ag qq — — 4CV(1 — x) with Cf = 4/3 and Agij = 
otherwise. This guarantees helicity conservation when 
the HVBM prescription for 75 in d dimensions is adopted 
to project onto definite helicity states |29| and defines 
the MS scheme in the polarized case [30(, which we use 
throughout our calculations. As a consequence of factor- 
ization, both the hard scattering cross sections and the 
parton distribution functions in Eq. ([7]) depend on the 
scale [if which is arbitrary, fij can be chosen differently 
than the renormalization scale fj, r at which ultraviolet sin- 
gularities are absorbed into the bare coupling and heavy 
quark mass. 

Note that dAa^ given in Eq. (Tiljj) receives an addi- 
tional singular contribution from the soft gluon parts of 
the diagonal splitting functions AP qq (x) and AP gg (x) 
in the factorization procedure, which is proportional to 
5(1 — x). Only then, all remaining singularities cancel in 

the sum of dAa^ and dAa^ , and the full expression for 
the subprocess cross section dAa a b in Eq. (|T2"j) is finite in 
the limit e — > 0. 

The numerical evaluation of J7]) for different measure- 
ment functions S can now be done in parallel with stan- 
dard Monte Carlo techniques by randomly generating a 
large sample of final-state configurations characterized by 
Xi,X2,x,y,6i, and Q 2 - The p and w-distributions regu- 
lating the singularities in dAa^' and the sum of dAa^ c± ^ 
and dAa c in Eqs. (TTS]) . tH}, and (TT6")) . respectively, need 
special attention. To this end, one inserts the definitions 
of distributions, given in Eq. (|30p of the Appendix, into 
Eqs. ([l"4"l) - (fl"l>]) and computes for each phase-space point 
a set of six correlated weights to account for all possi- 
ble configurations with x — 1 and y = ±1. The values 
of the measurement functions iS one is interested in are 
then multiplied by the appropriate weights and accumu- 
lated in different histograms. In principle, the choice 
for the parameters p £ [p, 1[ and w €]0, 2[ in (|30p does 
not matter as it only leads to different values for each of 
the individual, unphysical contributions at 0(a^) on the 
right-hand-side of Eq. (|12p but not for their sum. Large 
cancellations among the different terms in (TT2"|) can take 
place, however, if p is chosen too close to 1 or u> too close 

to [ali- 
To validate the numerical implementation of ([7]), we 
compare the results obtained with the Monte Carlo tech- 
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FIG. 1: Comparison of the results of our Monte Carlo 
code, dAaMC, with the analytical calculation dAa an of 
Refs. [Til l24l ] for single-inclusive charm production as a func- 
tion of transverse momentum and integrated over all ra- 
pidities. 



niques outlined above to those of the largely analyti- 
cal code developed and used in Refs. [l?], 0|- Since 
the analytical calculation is only applicable for single- 
inclusive heavy quark (or antiquark) hadroproduction, 
the comparison is done for charm production in longi- 
tudinally polarized pp collisions at a c.m.s. energy of 
= 200 GeV, without any experimental acceptance 
cuts. Figure [1] shows the difference of the numerical re- 
sults obtained with both codes, labelled as cIA<tmc and 
dAa an , as a function of the transverse momentum pj, of 
the charm quark, normalized to the analytical calcula- 
tion dAa an . As can be seen, deviations are at a level of 
a few per cent at most, well within the precision of the 
Monte Carlo integration for the relatively small sample 
of phase-space points used in Fig. [TJ 

Needless to say that all discussions in this Section also 
apply in the unpolarized case, and the corresponding 
expressions are obtained by appropriately replacing all 
helicity-dependent terms by their spin-averaged counter- 
parts. We fully agree with the results given in [23j |. 



III. PHENOMENOLOGICAL STUDIES 

A. Preliminaries 

Having laid out the technical framework in the previ- 
ous Section, we now turn to a detailed phenomenological 
study of heavy flavor hadroproduction in longitudinally 
polarized pp collisions and their semi-leptonic decays at 
RHIC. For comparison and to compute experimentally 



da 



(18) 



we also present results for the corresponding unpolarized 
quantities. We study the impact of the NLO QCD cor- 
rections on the polarized and unpolarized cross sections 
and quantify the theoretical uncertainties from different 
choices for unphysical factorization and renormalization 
scales, heavy quark masses, and parameters describing 
the hadronization of the heavy quarks. 

We concentrate on observables of immediate relevance 
for the RHIC spin program with collisions of longitu- 
dinally polarized protons at a c.m.s. energy of y/S = 
200 GeV. These are single-inclusive transverse momen- 
tum distributions of electrons and muons from semi- 
leptonic decays of charm and bottom quarks, and, in par- 
ticular, invariant mass spectra for two leptons observed in 
coincidence. Such measurements have been already car- 
ried out in spin-averaged pp collisions at RHIC [Til, EH 
and are intended with longitudinally polarized beams 
once sufficient statistics has been accumulated [l3j . 

We note that the leptons can stem from both charm 
and bottom quark decays which cannot be separated 
experimentally until displaced vertex detector upgrades 
have been installed. Therefore, our results always re- 
fer to the sum of charm and bottom production, their 
hadronization into D and B mesons, including c — > D, 
b — > B, and "cascade" b — > B — > D contributions, and 
the subsequent semi-leptonic decays of the heavy mesons 
into the observed leptons. We assume that electrons and 
muons are detected at central and forward rapidities, 
\rj e \ < 0.35 and 1.2 < < 2.2, respectively, which cor- 
responds to the acceptance of the PHENIX experiment 
Pi- 

The fragmentation of the heavy quarks into D and 
B mesons, i.e., D C ^ D and D b ^ B , is modeled by phe- 
nomenological functions extracted from fits to e + e~ data 
[H^ ] . D and B indicate a generic admixture of charm and 
bottom mesons. Contrary to fragmentation functions 
for light quarks and gluons into light mesons [l||, the 
non-perturbative functions describing the hadronization 
of heavy quarks are very hard, i.e., charm and bottom 
quarks only lose very little momentum when hadroniz- 
ing. The main effect of the fragmentation functions is to 
introduce a shift in the normalization of the heavy meson 
spectra. It depends mainly on the average momentum 
fraction z taken by the meson, while the details of the 
shape of D^^' Hq (z) have a negligible effect 19]. One 
can expect that ratios of cross sections, like in the exper- 
imentally most relevant double-spin asymmetry (|18[) . are 
much less affected by the actual choice of D®^ Hq (z). 
We use the functional form of Kartvelishvili-Likhoded- 
Petrov [20| with a single parameter ag controlling the 
hardness of 

D Q ^ H( i{z) = N Q z°"^(l-z), (19) 

where Nq = (olq + V)(olq + 2) to normalize the integral 
of dQ^ Hq (z) to one. We take a c = 5 and cvb = 15 from 
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Tab. 4 in Ref. as the default values in Eq. (fl9|) and 
vary them in the range 3 < a c < 7 and 10 < on, < 20, re- 
spectively, to estimate the uncertainties associated with 
the choice of «q . As in Ref. [lij] , the fragmentation is nu- 
merically performed by rescaling the heavy quark's three- 
momentum by z at a constant angle in the laboratory 
frame, i.e., ph q — zpq. The uncertainty introduced by 
this particular choice for the "scaling variable" z, which 
is not uniquely defined for D®~^ Hq , was shown to be not 
larger than scale and mass uncertainties [3l| and will be 
not considered further. 

The subsequent semi-leptonic decay of the D and B 
mesons into leptons is controlled by another set of phc- 
nomenological functions f H Q-+ e -f* which need to be ex- 
tracted from data as well. Here we use the spectra ob- 
tained in Ref. [3, H3] based on BaBar and CLEO data 
[HI ■ We note that we do not normalize our cross sections 
with the appropriate branching ratios for D — > e, B — > e, 
etc., which are all close to 10% J33J. Of course, branch- 
ing ratios drop out of experimentally relevant double-spin 
asymmetries (fig]) . 

The main motivation to study heavy flavor production 
with polarized beams at RHIC is the expected sensitiv- 
ity to the helicity-dependent gluon density through the 
tree-level gluon-gluon fusion process, gg — ► QQ, which is 
known to be dominant for unpolarizcd collisions up to the 
largest values of the heavy quark's transverse momentum 
currently accessible at RHIC [l5|, [r| . 

We will show, however, that the fractional contribu- 
tion of gluon-gluon fusion to the spin-dependent cross 
section depends crucially on the assumed set of polar- 
ized parton densities. Our default choice is the DSSV 
set Q, obtained in a global QCD analysis of the lat- 
est spin-dependent data, including those from RHIC on 
single- inclusive pion and jet production 0, 0]. Due to 
the smallness of Ag(x) in the DSSV set and a node in 
the x-shape near x ~ 0.1 [8|, the qq annihilation sub- 
process, qq — > QQ, turns out to be the dominant mech- 
anism for charm and bottom production in polarized pp 
collisions at RHIC. This is in sharp contrast to naive 
expectations based on unpolarized results. For compari- 
son and to study the sensitivity to Ag(x), we adopt also 
two alternative sets of spin-dependent parton densities, 
GRSV(std) [H and DNS(KRE) [H, both characterized 
by a positive gluon polarization of moderate size. In 
general, for gluon polarizations from current QCD fits 
[3; HH, [35| , the double-spin asymmetries for leptons from 
heavy flavor decays all turn out to be small, often well 
below the one percent level, making their measurement 
very challenging. This is in particular true for single- 
inclusive lepton observables; see below. 

Heavy flavor production at RHIC cannot compete with 
the statistical precision achievable for more abundant 
probes of the nucleon's spin structure, like pions and jets 
[5|, Q , which are already used in global fits f8] . Neverthe- 
less, measurements of double-spin asymmetries related to 
heavy flavor production will be crucial for further test- 
ing and establishing the all important concept of factor- 



ization and universality for helicity-dependent scattering 
processes and parton densities, respectively. The under- 
lying dynamics of the partonic scattering processes, i.e., 
gg — > QQ and qq — > QQ, is very much different as com- 
pared to the multitude of QCD processes driving the pro- 
duction of light hadrons [3(| or jets [13] ■ 

In the computation of the NLO unpolarized cross sec- 
tions in (|18[) . which proceeds along similar lines as out- 
lined in Sec. II, for details, see Ref. (23|, we use the NLO 
CTEQ6M parton densities [38[ and values for the strong 
coupling a s . Since the DSSV analysis [|[ does not provide 
a LO set of spin-dependent parton distributions, our LO 
results always refer to the Born part of the full NLO cal- 
culation, i.e., they are computed with NLO parton densi- 
ties and values for a s . Strictly speaking this is, of course, 
inconsistent as it introduces some unwanted scheme de- 
pendence into a tree-level quantity. Nevertheless, the LO 
results should give a faithful estimate of the relevance of 
NLO corrections. As will be demonstrated below, they 
turn out to be sizable and rather different for unpolarized 
and polarized cross sections such that they do not cancel 
in experimentally relevant double-spin asymmetries. 

We take m c = 1.35 GeV and mj, = 4.75 GeV as ref- 
erence values for the charm and bottom quark mass 
and vary them in the range 1.2 < m c < 1.5 GeV and 
4.5 < nib < 5.0 GeV, respectively, to estimate the result- 
ing mass uncertainties. For the factorization and renor- 
malization scale we take £*/ = fi r = £( m Q + [(Pt) 2 + 

(p^) 2 ]/2) 1 / 2 with £ = 1 as the central value. As usual, 
the sensitivity of the cross section to missing higher order 
corrections is estimated by varying fj, f and fj, r simultane- 
ously in the range 1/2 < £ < 2. Following the procedure 
used for unpolarized charm and bottom production at 
RHIC in Ref. [l4[ , we also vary fi f and fi r independently 
in the same range of £ and combine the ensuing uncer- 
tainty with the one stemming from variations of m c j, in 
quadrature. Unless stated otherwise, we use the central 
values for /i/, fi r , m Ct b, and a c ^ given above. 

B. Heavy Flavor Cross Sections And Correlations 

We begin our detailed numerical studies with a dis- 
cussion of unpolarized and polarized cross sections for 
various decay lepton distributions accessible at RHIC. 

Figure [5] shows the single-inclusive transverse momen- 
tum spectrum of electrons from charm and bottom de- 
cays in LO and NLO accuracy, integrated over the angu- 
lar acceptance of the PHENIX detector, i.e., \r/ e \ < 0.35 
PH . Similar results are obtained for the STAR experi- 
ment [l2j with its larger acceptance for electrons at cen- 
tral rapidities, \r] e \ < 1, and hence not shown. The 
transverse momentum p^, is limited to a region which 
should be accessible with luminosities envisaged in lon- 
gitudinally polarized pp collisions at RHIC. Photon con- 
version, 7 — > e + e~, and ir° — > je + e~ Dalitz decays are 
the dominant source of electron background for such mea- 
surements and may require an additional cut p^ > 1 GeV 
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FIG. 2: Scale dependence of the single-inclusive transverse 
momentum spectrum of electrons from heavy quark decays 
at central rapidities \r) e \ < 0.35 in unpolarized (upper panel) 
and polarized (middle panel) pp collisions at RHIC. All scales 
are varied simultaneously, i.e., fXf = [i r = £(w.q + [(p^) 2 + 
(p!p) 2 ]/2) 1/2 , in the range 1/2 < £ < 2 (shaded bands). The 
solid lines correspond to the default choice £ — 1. Note that 
the LO results are rescaled by a factor of 0.1, and all polarized 
cross sections are multiplied by —1. The lower panel shows the 
ratio of NLO to LO polarized and unpolarized cross sections 
(A"-factor). 



[ill H2j . Recall that the branching ratios of about 10% 
are not included in the cross sections shown in Fig. [2j 

The solid lines are obtained with the default values of 
the heavy quark masses, scales, parameters, and parton 
densities stated in the previous Subsection. The shaded 
bands indicate the theoretical uncertainty from varying 
the factorization and renormalization scale simultane- 
ously in the range 1/2 < £ < 2 specified above. Note 
that the LO results are rescaled by a factor of 0.1, and 
all polarized cross sections are multiplied by -1 to dis- 
play them on a logarithmic scale. The bottom panel of 
Fig. [5] gives the resulting unpolarized and polarized "K- 
factors" , defined as usual by the ratio 



^ _ d(A)a NLO 
K — 



d{A)a 



LO 



(20) 



One notices that the NLO corrections are sizable in the 
unpolarized case, K ~ 3, but moderate for polarized pp 
collisions, except for the region p^, > 5 GeV. Here, the 
polarized cross section approaches a node, and perturba- 
tive corrections are artificially enhanced. 

Less pronounced NLO corrections for polarized cross 



FIG. 3: Same as in Fig. [2] but for the single-inclusive trans- 
verse momentum spectrum of muons from heavy quark decays 
at forward rapidities 1.2 < |^ M | < 2.2. 



sections are a rather generic feature and have been ob- 
served already for other hadronic processes such as single- 
inclusive pion [36| and jet [37j production. To some ex- 
tent this behavior can be traced back to the less singular- 
scale evolution of polarized parton densities at small mo- 
mentum fractions x [30| . This has the effect that the 
partonic threshold region, which is the source of large 
logarithmic corrections associated with the emission of 
soft gluons, is less emphasized in the convolution (J7J than 
in the unpolarized case. Specifically for heavy flavor pro- 
duction, it was noticed in [l7| that large NLO corrections 
to the gluon-gluon fusion process related to amplitudes 
with a gluon exchange in the i-channel are independent 
of the helicities of the interacting gluons and hence do 
not contribute to the polarized cross section. Substan- 
tially different i^-factors for unpolarized and polarized 
cross sections immediately imply that Born level esti- 
mates for double-spin asymmetries (|18[) can serve only 
as very rough estimates. In general, they are insufficient 
for any quantitative analysis such as a global QCD ex- 
traction of spin-dependent parton densities. 

As can be inferred from Fig. [21 theoretical ambiguities 
due to the choice of [i r and /jj in J7| are reduced in the 
polarized but not in the unpolarized case. Along with 
the observation of large QCD corrections, this indicates 
the need for next-to-next-to-leading order corrections for 
the unpolarized cross section to better control the depen- 
dence on fj,f and fx r . We note that substantial progress 
toward this direction has already been made recently 
[39l | . mainly to allow for precision studies with the large 
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FIG. 4: Same as in Fig. [5] but for the invariant mass spec- 
trum of electron- muon correlations from heavy quark decays. 
Electrons are detected at central rapidities, \r) e \ < 0.35, and 
muons are detected at forward rapidities, 1.2 < < 2.2. 
The transverse momenta of both electrons and muons are re- 
quired to be larger than 1 GeV. 



FIG. 5: Same as in Fig.[2]but for the invariant mass spectrum 
of muon- muon correlations from heavy quark decays. Both 
muons are detected at forward rapidities, 1.2 < |tj m | < 2.2, 
but in different hemispheres. The transverse momentum of 
each muons is required to be larger than 1 GeV. 



amount of top quarks expected to be produced at the 
CERN-LHC. In addition, fixed order calculations need 
to be amended by all-order resummations if In p®/m,Q 
becomes large. This was achieved in 2l| but is not re- 
ally relevant for our discussions here since we are mainly 
interested in the region where p® ~ ttlq. We postpone a 
discussion of theoretical uncertainties due to the choice of 
rriQ and o>q in Eq. (TIT))) , as well as the effect of varying fi r 
and \ij independently, until the end of this Subsection. 

The single- inclusive transverse momentum spectrum of 
muons from heavy quark decays is shown in Fig. [3J in LO 
and NLO accuracy. The pseudorapidity 77^ of the muon is 
integrated in the range 1.2 < \r\p\ < 2.2 corresponding to 
the angular acceptance of the PHENIX experiment. All 
observations made in Fig. [2] regarding the relevance of 
NLO corrections, the behavior of the if-factor, and the 
dependence on fj,f >r apply also here. The polarized K- 
factor stays even closer to one than in Fig. [2] as dAcr/dpt^ 
develops no node in the p 1 ^ range shown. 

By glancing at the relative sizes of the unpolarized and 
polarized single-inclusive transverse momentum spectra 
for electrons and muons shown in Fig. [5] and [3J re- 
spectively, it becomes immediately obvious that the cor- 
responding double-spin asymmetries A e LL and A^ L , to 
which we turn to in Subsec. IIIIDl are very small if the 
most up-to-date DSSV parton densities [|[ are used. Ex- 
pected asymmetries of the order of a few tenths of a 



percent are extremely challenging experimentally as sys- 
tematic uncertainties, like from the determination of the 
relative beam luminosities at RHIC, are of similar size 
Q. At higher p^, where double-spin asymmetries are 
largest, the single-inclusive cross sections in Figs. [5] and 
[3Jhave dropped already several orders of magnitude from 
their peak values, and measurements require substantial 
integrated luminosities. 

More promising appear to be observables where both 
the heavy quark and the heavy antiquark decay semi- 
leptonically, and both leptons are observed in coinci- 
dence. This is also where our numerical phase-space in- 
tegration and the flexible Monte Carlo code for polarized 
heavy flavor hadroproduction introduced in Sec. [TT] be- 
come truly essential. Particle correlations are hard, and 
often impossible, to compute at NLO with largely ana- 
lytical methods, see, e.g., [40| . 

Figure U shows our results for the invariant mass spec- 
trum of electron-muon correlations from semi-leptonic 
decays of D and B mesons within the angular accep- 
tance of the PHENIX detector, i.e., |?y e | < 0.35 and 
1-2 < \t)h\ < 2.2. In addition, we require a minimum 
transverse momentum for both leptons of 1 GeV as re- 
quired by experiment. As in Figs. [5] and [3J results shown 
as solid lines are obtained with the default choice of pa- 
rameters. Again, shaded bands indicate the theoretical 
uncertainty from varying the factorization and renormal- 
ization scales simultaneously in the range 1/2 < £ < 2. 
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FIG. 6: Scale and mass uncertainties for the polarized in- 
variant mass spectrum of electron-muon correlations at NLO 
accuracy, using the same cuts as in Fig. [3] The upper and 
lower panel shows the contribution from cc and bb decays, re- 
spectively, multiplied by —1. The dashed lines indicate the 
range of uncertainties for /i/ = /j, r , 1/2 < £ < 2, and central 
values of m Cj i, as in Fig. [4] The solid curves are for £ = 1. 
The effect of varying fif and /i r independently in the same 
range of £ for fixed ra c = 1.35 GeV, mi, = 4.75 GeV and for 
1.35 < m c < 1.5 GeV, 4.5 < m b < 5.0 GeV is illustrated by 
hatched and solid bands, respectively. In the latter case, scale 
and mass uncertainties are combined in quadrature (see text). 



The if-factor, shown in the lower panel of Fig. 0] is 
smaller than what was found for the single-inclusive ob- 
servables in Figs. [2] and in the unpolarized case. Still, 
NLO corrections differ considerably for the unpolarized 
and polarized invariant mass spectra. Again, the correc- 
tions are such that the corresponding double-spin asym- 
metry is reduced at NLO accuracy. The scale uncertainty 
is significantly smaller for the spin-dependent cross sec- 
tion with NLO corrections included. The improvement 
in the helicity-averaged case is much less pronounced. 

Compared to the single-inclusive results in Figs. [2] 
and [31 the cross sections obtained for the electron-muon 
invariant mass spectrum are smaller, but d(A)a / dm e ^ 
drops much less with increasing m eAI than d(A)a / dp^ 
with increasing pfp or pt^. This makes measurements of 
A e £ L at comparatively large values of m efi feasible. We do 
not consider here correlations with back-to-back electrons 
at central rapidities. Electron-muon correlations are phe- 
nomenologically more interesting due to their asymmet- 
ric kinematics with respect to rapidity, probing the in- 
teracting partons at different momentum fractions x as 
will be demonstrated below. 

The corresponding invariant mass spectrum for two 
muons from D and B meson decays observed in coin- 
cidence is shown in Fig. Both muons are required to 
have 1.2 < \n^\ < 2.2 and pf^ > 1 GeV, with one muon 



detected at forward (positive) and one muon detected at 
backward (negative) pseudorapidities. Again, the cross 
sections d(A)a decrease rather slowly with increasing in- 
variant mass m MM . This observable is very demanding in 
terms of required Monte Carlo statistics as can be seen 
by the still fairly pronounced fluctuations, most notice- 
able in the unpolarized if- factor. The general trend and 
features of the cross sections are, however, reliable. As 
before, NLO corrections are more significant in the unpo- 
larized case, where K ~ 2. Unfortunately, the reduction 
of the theoretical ambiguities related to the choice of fif, r 
is only marginal at NLO. 

We now turn to a more detailed discussion of theoreti- 
cal uncertainties for the observables discussed in this Sec- 
tion, taking the phenomenologically interesting invariant 
mass spectrum for electron-muon correlations, presented 
in Fig. 21 as an example. Qualitatively very similar re- 
sults are obtained for the other cross sections given in 
Figs. [2] [3l and [5] and hence not shown here. 

The impact of varying p,f and p, r independently 
is shown in Fig. [6] for the electron-muon invariant 
mass spectrum in polarized pp collisions at NLO 
accuracy. Since we are also interested in variations 
of m c ,b, the contribution from cc and bb decays are 
shown in separate panels and add up to —dAa/dm etl 
discussed in Fig. Ill Following Ref. jlij], we com- 
pute our results for seven different settings of scales 

M/,r = J/„K + [(P?) 2 + b?) 2 ]/2) 1 / 2 , nsing = 
{(1, 1), (2, 2), (1/2, 1/2), (1, 1/2), (2, 1), (1/2, 1), (1, 2)} 
and keeping tuq fixed to their central values 
m c = 1.35 GeV and nib — 4.75 GeV. The envelope 
of all resulting curves defines the scale uncertainty and 
is show as hatched bands in Fig. [6] For comparison, 
the dashed lines indicate the range of uncertainties for 
standard choice fj,f — fi r used in Fig.[U As can be seen, 
taking p,f =/= p r does not significantly enlarge the scale 
ambiguities for the polarized cross section, in particular, 
for the contribution from bottom quarks shown in the 
lower panel. 

The solid bands in Fig. [5] take also variations of 
m c and nib into account. The recipe we follow here 
is again similar to the one used to estimate theoreti- 
cal uncertainties for unpolarized heavy flavor produc- 
tion [3, [3l|. In practice, we add scale and mass un- 
certainties in quadrature, and the envelope of all re- 
sults is defined by C + [(M+ - C) 2 + {M+ Q - C) 2 ] 1/2 
and C -\{C - M-f + (C - M mQ ) 2 ] 1/2 - Here, C de- 
notes the results obtained for central values of scales and 
masses. (M~) are the maximum (minimum) cross 

sections computed for p,f ^ /i,., m c — 1.35 GeV, and 
nib = 4.75 GeV, as depicted by the hatched bands. Cor- 
respondingly, M+ Q (M~ Q ) denote the maximum (min- 
imum) cross sections for £f^ r = 1 and varying m c and 
rrifc in the range 1.35 < m c < 1.5 GeV and 4.5 < m& < 
5.0 GeV, respectively. In general, the combined uncer- 
tainties are much smaller for bb than for cc production 
and decays, which is not too surprising. In both cases, 
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C. Subprocess, Charm, and Bottom Fractions 
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FIG. 7: Dependence of the polarized (solid lines) and unpolar- 
ized (dashed lines) invariant mass spectra for electron-muon 
correlations at NLO accuracy on the choice of fragmentation 
parameters a c (upper panel) and ab (lower panel) defined in 
Eq. (|19[) . Displayed are the relative deviations for cc and 
bb decays using a c = 3, 7 and otb = 10, 20 with respect to 
the cross sections obtained for our default values a c — 5 and 
at = 15 respectively, [if, fj, r , and heavy quark masses 
are taken at their central values, and the same experimental 
cuts as in Fig. [3] are adopted. 
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variations of toq add noticeably to the theoretical uncer- 
tainties. 

The dependence of the cross sections on the choice of 
cvq in the non-perturbative function dQ~* Hq (z) describ- 
ing the hadronization of the heavy quarks into D and 
B mesons, see Eq. (fT9|) . is illustrated in Fig. [7] Again, 
we take the invariant mass spectrum for electron-muon 
correlations as an representative example. We vary a c 
and ab in the range [19| 3 < a c < 7 and 10 < a& < 20, 
respectively, and show the impact on the invariant mass 
spectrum as relative uncertainty with respect to the re- 
sults obtained for the central values a c = 5 and ab = 15 
used in Fig. [4] 

It turns out that polarized and unpolarized invariant 
mass spectra arc affected very much in the same way 
by variations of ol c ^- For charm production and tak- 
ing 3 < a c < 7, it roughly amounts to a shift in the 
normalization of the cross sections by ± 20 4- 30%. The 
impact of varying ab on the contribution to d(A)<r / dm e ^ 
from bottom decays is significantly smaller, up to about 
± 10% deviation from the results for ab = 15, but is less 
uniform with m ell . These observations have the impor- 
tant implication that theoretical uncertainties associated 
with the actual choice of a c ^ drop out to a large extent 
for experimentally relevant double-spin asymmetries All 
discussed in Sec. MIDI below. 



FIG. 8: Fractional amount of different partonic subprocesses 
at NLO accuracy (left column) and of charm, bottom, and 
cascade (b — > c) decays (right column) contributing to the 
single-inclusive transverse momentum spectrum of electrons 
shown in Fig. [2] Results are shown for unpolarized (upper 
row) and polarized (middle and lower rows) pp collisions at 
RHIC using the CTEQ6 [H, DSSV [|, and GRSV [3 set 
of parton densities, respectively. 

We now take a detailed look at the fractional contribu- 
tions of the different partonic hard scattering processes 
to the cross sections shown in Figs. [5]- [5] This will help 
to understand the dependence of the double-spin asym- 
metries on different sets of polarized parton densities, to 
be discussed in the next Subsection. Since charm and 
bottom decays both contribute to the lepton spectra, we 
also present their fractional contributions. This includes 
also the "cascade" decay b — * c — * e, fi, which is modeled 
following the procedure discussed in [14]]. It is found to 
be negligible for all observables we are interested in. 

The left-hand-side of Fig. [5] shows the contributions of 
the three different subprocesses at NLO accuracy, with 
gg, qq, and qg initial-states, to the single-inclusive decay 
electron spectra shown in Fig. [2] In the unpolarized case 
(upper panel) , gluon-gluon fusion is the by far dominant 
subprocess for heavy flavor production at RHIC energies, 
with qq annihilation becoming somewhat more relevant 



12 



at larger values of transverse momentum p T . Interest- 
ingly enough, the genuine NLO, i.e., a s suppressed, qg 
scattering process also contributes very significantly at 
larger p T , even exceeding the qq annihilation cross sec- 
tion. This observation can be linked to the abundance 
of gluons at all momentum fractions x [381 ] . This im- 
plies that the partonic flux relevant for qg scattering, i.e., 
q(xi, fif)g(x2, [if), is much larger than the corresponding 
flux for qq annihilation, in particular, at the medium-to- 
large momentum fractions X\^. relevant for RHIC. This 
compensates for the 0{a s ) suppression in the qg hard 
scattering channel. In pp scattering, e.g., at the TeVa- 
tron, where antiquaries are "valence" quarks in the an- 
tiproton beam, this is different, and the qq flux is much 
enhanced. A similar observation concerning the relevance 
of the qq annihilation channel was made also for fixed- 
target experiments in Ref. [24j . where it is expected to 
contribute very significantly to charm hadroproduction 
in proposed pp collisions at the GSI-FAIR facility, but 
not in pp scattering of similar c.m.s. energy planned at 
J-PARC. 

In general, the situation is much more involved in the 
polarized case, where both hard scattering cross sections 
and parton densities are not positive definite and can 
contribute with either sign, depending on the kinematics 
relevant for a particular process. In the vicinity of sign 
changes, large cancellations are to be expected. As we 
shall demonstrate below, depending on the chosen set of 
polarized parton densities, the subprocess fractions can 
differ considerably from each other and often gluon-gluon 
fusion does not dominate, in contrast to the unpolarizcd 
case. 

The middle panel of Fig. [5] shows our results for the po- 
larized subprocess fractions obtained with the DSSV set 
Q, our default choice of parton densities used in Figs.O- 
[5J At small p T , the cross section is entirely dominated 
by qq annihilation, contrary to the unpolarized case. To- 
wards larger p T , both gg and qg processes contribute sig- 
nificantly but with opposite sign, leading to strong can- 
cellations. This happens, however, in a kinematic region 
close to a sign change of the cross section at p T ~ 7 GeV. 

A rather different pattern of fractional subprocess con- 
tributions can be found in the bottom panel of Fig. [SJ 
where the GRSV(std) parton densities were used. 
Note that the cross section has a sign change near p T = 
2 GeV. This explains the complicated behavior of the ra- 
tios in this region and makes it very awkward to display 
them properly. Like in the unpolarized case, gluon-gluon 
fusion is the most important contribution to the cross 
section. At larger p T , qq and qg subprocesses become 
more relevant, both contributing with the opposite sign 
than gg scattering, leading again to fairly significant can- 
cellations. 

The gross features of the results in Fig. [SJ obtained 
with DSSV and GRSV parton densities can be readily 
understood by comparing the size and sign of the in- 
dividual quark, antiquark, and gluon densities in both 
sets, see, e.g., Fig. 2 in Ref. [8j. Since the decay electrons 
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FIG. 9: Same as in Fig. [8] but now for the single-inclusive 
transverse momentum spectrum of muons shown in Fig. [3] 

stem from heavy (anti)quarks produced at central rapidi- 
ties, the interacting partons have very similar momentum 
fractions, i.e., x\ ~ X2- Therefore, Ag(x±)Ag(x2) > 0, 
irrespective of the node in the DSSV gluon distribution, 
and the sign of the gg contribution follows the sign of the 
hard scattering cross section, which changes from posi- 
tive at small px to negative at larger px values. 

Due to helicity conservation, dAa q q = —da q q < 
[T3|, and the sign of the qq contribution depends on 
the individual parton densities for each quark and an- 
tiquark flavor. All sets of polarized parton densities have 
Au(x) > and Ad(x) < 0, resembling the features of 
the naive quark model. The GRSV(std) set [34[ assumes 
an SU (3) symmetric sea, with all antiquark polarizations 
being negative, such that for the dominant u-quarks one 
has Au(xi)Au(x2) < 0, resulting in a net positive con- 
tribution to the cross section. This is exactly opposite 
in the DSSV set [8J, where Au(xi)Au(x2) > 0, unless X2 
gets very large. The genuine NLO qg subprocess cross 
section, as well as the sum of all quark and antiquark 
polarizations, AE = ^2 q [Aq + Aq], are both positive. 
This implies that the sign of the qg contribution de- 
pends on the sign of Ag(x) in the relevant region of x, 
which turns out to be positive for both GRSV(std) and 
DSSV. Overall, the fractional contributions of the indi- 
vidual subprocesses to the single-inclusive decay electron 
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FIG. 10: Same as in Fig. [8] but now for the invariant mass 
spectrum of electron-muon correlations shown in Fig. [4] 



FIG. 11: Same as in Fig. [10] but now for the invariant mass 
spectrum of muon-muon correlations shown in Fig. [5] 



spectrum are essentially controlled by the modulus of the 
polarized gluon density, |Ag(:r)|, which is much larger for 
the GRSV(std) set, i.e., \Ag{x)\ GRS v > \Ag{x)\ DSS v- 
The bigger the gluon density, the closer the result is to 
what we have found in the unpolarized case. It turns 
out that even for the moderate gluon polarization of the 
GRSV(std) set, the gluon-gluon channel prevails for all 
Pt values shown in Fig. [H] 

The right-hand-side of Fig. [5] shows the fractional con- 
tributions of the charm, bottom, and "cascade" b — > c 
decays to the single- inclusive transverse momentum spec- 
trum of electrons. For pf, < 2GeV, almost all electrons 
originate from charm decays, but above the bottom con- 
tribution catches up, yielding about 25% at pj* = 6 GeV 
in the unpolarized case shown in the upper panel of 
Fig. [8] Eventually, at somewhat larger values of pfi, it 
becomes dominant, as was shown in [14] . 

As for the subprocess fractions, the corresponding re- 
sults in the polarized case depend on the choice of parton 
densities, and cancellations and possible sign changes of 
the individual hard scattering subprocesses dAa a b fur- 
ther complicate their interpretation. The only feature 
common to both the unpolarized and the polarized in- 
clusive electron spectra is the smallness of the b — > c — ► e 
cascade decay contribution. The results obtained with 
the DSSV set (middle panel) show a very sizable bot- 



tom contribution, exceeding the c — > e decay above 
Pj, ~ 4 GeV. However, this is due to a sign change of 
the c — > e cross section at p^, ~ 7 GeV, and above, c — > e 
and b — > e contribute on equal footing. Choosing the 
GRSV(std) distributions instead (lower panel), both the 
c — * e and the b —> e cross sections change from positive 
to negative at 2 and 4 GeV, respectively, with b — > e start- 
ing to be the dominant contribution above p^ ~ 6 GeV. 

The fractional contributions of the different partonic 
hard scattering processes (left-hand-side) and heavy fla- 
vor decays (right-hand-side) to the single-inclusive muon 
cross section shown in Figs. [3] are given in Fig. [9l The 
results are qualitatively very similar to the ones depicted 
in Fig. [5] and discussed above. Again, the polarized 
subprocess fractions very much depend on the choice of 
parton densities, and the interpretation is obscured by 
sign changes and large cancellations among the differ- 
ent contributions. Compared to the single-inclusive elec- 
tron spectrum at central rapidities, bottom decays con- 
tribute less to the muon transverse momentum spectrum 
at 1.2 < | r^l < 2.2, even up to Pj, — 10 GeV. Its contri- 
bution is rather flat with respect to p^ and amounts to 
about 40% (20%) in the (un)polarized case. 

Corresponding results for the electron-muon and 
muon-muon invariant mass distributions, shown in Fig.0] 
and O can be found in Fig. [TU] and [TT1 respectively. As 
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FIG. 12: Double-spin asymmetry A LL for single-inclusive 
electrons from charm and bottom decays at RHIC, computed 
at NLO accuracy for three different sets of polarized par- 
ton densities: DSSV [| (solid line), GRSV(std) (dashed 
line), and DNS (KRE) [3l] (dot-dashed line). Electrons are 
restricted to central rapidities \r) e \ < 0.35. 



before, subprocess fractions can be found on the left- 
hand-side and contributions from different heavy flavor 
decays on the right-hand-side of the plots. Note that in 
all panels of Figs. [10] and [Til the bin corresponding to 
the smallest invariant mass only has a small number of 
entries due to the cuts > 1 GeV, and the numeri- 
cal results for that bin should be taken with caution. In 
general, cancellations among different subprocesses are 
found to be less pronounced in Figs. fTUl and [TTl except 
for small invariant masses, say, below 4 GeV, where a sign 
change in the polarized cross section occurs. 

Gluon-gluon fusion is even more dominant for electron- 
muon and muon-muon correlations than for single- 
inclusive decay lepton observables, with both qq and qg 
subprocesses being negligible in the unpolarized case (up- 
per row). This is also the case for the polarized cross 
section if the GRSV (std) parton distributions are cho- 
sen. For the DSSV set, qq annihilation remains domi- 
nant for electron-muon correlations, but gluon-gluon fu- 
sion contributes significantly to muon-muon correlations 
for to w > 10 GeV. As will be shown in the next 
Subsection, back-to-back muon-muon correlations with 
1-2 < 1 77^ | < 2.2 probe on average fairly large momentum 
fractions, (x) > 0.1, where the DSSV Ag(x) is positive 
and larger than the sea quark polarizations, such that 



FIG. 13: Same as in Fig. [12] but for the single-inclusive muon 
spectrum at forward rapidities 1.2 < \r]^\ < 2.2. 



dA<7gg > dAa q q. Since the DSSV Au(x) turns negative 
at large x, there are additional cancellations among the 
different quark flavors in the qq annihilation channel, as 
Au(xi)Au(x 2 ) < and Ad(x 1 )Ad(x 2 ) > 0. 

As in Figs. [8] and [9] the contribution from the cascade 
b — > c — > e,fi decay is found to be negligible for lepton- 
lepton correlations. In particular, the muons in back-to- 
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FIG. 14: Same as in Fig. 1121 but for the invariant mass spec- 
trum for electron-muon correlations. Electrons and muons are 
restricted to |?7e| < 0.35 and 1.2 < < 2.2, respectively. In 
addition, we demand Pt > 1 GeV. 



back correlations originate mainly from charm decays, as 
can be seen on the right-hand-side of Fig.[TTJ At the same 
invariant lepton- lepton mass, bottom quark decays con- 
tribute more significantly to electron-muon correlations 
than to muon-muon correlations for m eM , m w > 5 GeV. 



D. Double-Spin Asymmetries 

The quantities of actual interest in experiments ex- 
ploiting polarized beams and targets are the double-spin 
asymmetries All defined in Eq. (|18[) . Experimental nor- 
malization uncertainties conveniently cancel to a large 
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FIG. 15: Same as in Fig. [12] but for the invariant mass spec- 
trum for muon-muon correlations. The muons are restricted 
to forward rapidities 1.2 < |7j M | < 2.2 and have to be in op- 
posite hemispheres. In addition, we demand > 1 GeV. 

extent in the ratio (fT5)) . In general, this does not happen 
for higher order QCD corrections or the various sources 
of theoretical ambiguities as we shall demonstrate below. 
Nevertheless, it is often erroneously assumed that LO es- 
timates for All give reliable results which can be used 
in quantitative QCD analyses. 

We start by giving expectations for various double-spin 
asymmetries at NLO accuracy in Figs. [T^]- [HI based on 
the polarized and unpolarized cross sections for decay 
lcpton transverse momentum and invariant mass distri- 
butions presented in Figs. [2]- Apart from our default 
choice of DSSV polarized parton densities [8| , which leads 
to very small asymmetries throughout, we adopt also two 
alternative sets, GRSV(std) |3 and DNS(KRE) 
Unlike DSSV, both sets are characterized by a positive 
gluon polarization of moderate size and an almost SU (3) 
symmetric sea. We refrain from using outdated models 
with a large, but strongly disfavored gluon polarization in 
the x range already probed by RHIC pp and fixed target 

data Ha MUg. 

With the exception of the double-spin asymmetry A LL 
for single-inclusive electrons from charm and bottom de- 
cays shown in Fig. [TSJ differences in the results obtained 
with GRSV(std) and DNS(KRE) parton distributions 
are readily explained by the slightly larger Ag(x) in the 
GRSV set. For A e LL , the result based on the DNS(KRE) 
set is strongly affected by cancellations between the gg 
subprocess on the one hand, and the qq, qg processes on 
the other hand, leading to an essentially zero spin asym- 
metry in the pj, range shown. Cancellations among the 
different subprocesses contributing with different sign are 
less pronounced for the GRSV set, see the lower left panel 
of Fig. [51 due to a significantly less negative Au density 
at x ~ 0.1. Cancellations in conjunction with the small- 



ness of A LL for all p? make this observable not really 
suited for studies of the nucleon's spin structure. 

The double-spin asymmetry A^ LL for single-inclusive 
decay muons at 1.2 < | r/^ | < 2.2, presented in Fig. [T3l 
shows a much better correlation of the size of Ag(x) 
and A^ L . As was demonstrated in Fig. [5] cancellations 
among the different subprocesses are less pronounced 
than for A e LL . For the same value of transverse momen- 
tum, the obtained spin asymmetries are about a factor 
of two larger for A^ L (pf^,) than for A LL (j3j,). However, 
they are still significantly smaller than projections based 
on very large positive, but outdated gluon polarizations 
like GRSV(max) H, see, e.g., Fig. 3 in Ref. 0. Never- 
theless, with sufficient statistics accumulated, spin asym- 
metries of 0(0.5%) should be measurable. Based on the 
cross sections given in Figs. [2] and [3] we estimate that an 
integrated luminosity of about 1 fb _1 is required. 

The best suited observables related to heavy flavor 
hadroproduction in longitudinally polarized pp collisions 
at RHIC are double-spin asymmetries for invariant mass 
spectra of electron-muon and muon-muon correlations 
shown in Fig. [14] and 1151 respectively. At small invari- 
ant mass, the corresponding cross sections are smaller 
than for single-inclusive transverse momentum distribu- 
tions at similar values of but fall off much slower 
with increasing invariant mass, cf. Figs. [2] - [5] This 
should allow for measurements of AV\ and A'i'i up to 
TOe^i.MM = 10 -r- 15 GeV, where spin asymmetries can be 
up to the I-t-2 % level for the GRSV(std) and DNS(KRE) 
parton densities. An integrated luminosity of a few hun- 
dred pb _1 should be sufficient. In addition, the size 
of A^j^ 1 ^ and A<?(x) in the relevant range of momen- 
tum fractions x, see below, are nicely correlated. In the 
absence of shifted vertex detectors at the RHIC experi- 
ments, A^t is the observable with the cleanest sample of 
charm decays for all m MM shown in Fig. [TS] and irrespec- 
tive of the set of polarized parton densities used in the 
calculations, cf. Fig. [11] 

Figure clearly illustrates the inadequacy of com- 
puting double-spin asymmetries based on LO estimates 
of heavy flavor cross sections. Depicted is the i^-factor, 
i.e., All(NLO)/All(LO), for all double-spin asymme- 
tries presented in Figs. [T2] - EH computed in each case 
using the DSSV polarized and CTEQ6 unpolarized par- 
ton densities. Other sets of polarized parton densities 
yield qualitatively very similar results. On average, LO 
estimates for All are about a factor of two larger than 
corresponding calculations at NLO accuracy and depend, 
in case of the single- inclusive observables A LL and A^ L , 
strongly on p e j^ . This reflects the difference of -FGfactors 
for the polarized and unpolarized cross sections found in 
Figs. [2] - [5] and invalidates any approximation based on 
constant iGfactors or the idea that higher order QCD 
corrections cancel in All- 

As was already illustrated in Fig. [7] theoretical uncer- 
tainties associated with the actual choice of the param- 
eters a Ct b in the non-perturbative function D®^ Hq (z) 
given in Eq. (fTT)]) cancel to a large extent in double- 
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FIG. 16: Impact of NLO QCD corrections on the double-spin 
asymmetries shown in Figs. [p2]- US] Depicted is the K- factor, 
i.e., All(NLO)/All(LO), computed in each case using the 
DSSV polarized and CTEQ6 unpolarized parton densities. 



spin asymmetries. Unfortunately, this is not the case 
for ambiguities related to the choice of scales /i/. r - As 
we have discussed in Sec. IIIIB1 the dependence of un- 
polarized heavy flavor cross sections on variations of 
fj,f >r is in general more pronounced than in the polar- 
ized case, see, e.g., Fig. [4] for electron- muon correlations. 
This can cause sizable ambiguities also for ratios of cross 
sections, like double-spin asymmetries. As a represen- 
tative example, we show in Fig. [T7] the dependence of 
the double-spin asymmetry A e ^ L for electron-muon cor- 
relations on variations of /x/ ir . The shaded band in the 
upper panel of Fig. [T71 illustrates the uncertainty on A e ^ L 

if Hf = n r = £(mjg + [{p®) 2 + (p^) 2 ]/2) 1 / 2 are varied si- 
multaneously in the range 1/2 < £ < 2. The lower panel 
gives the relative deviation of A e ^ L for £ = 1/2, 2 with 
respect to the result obtained for our default value £ = 1. 

The scale uncertainties are quite substantial and not 
uniform as a function of the invariant mass m e(1 . Never- 
theless, the asymmetries obtained with the DSSV parton 
densities are still much smaller than for sets with larger 
gluon polarizations, like DNS(KRE) or GRSV(std), as 
can be inferred by comparing with the results given in 
Fig. [14] Qualitatively similar effects as in Fig. [17] are 
found for the other double-spin asymmetries discussed in 
this Subsection. We refrain from varying /i f and fi r inde- 
pendently, which increases the uncertainties only slightly 
for the specific observables and kinematics we are inter- 
ested in, cf. Fig. (5] Similar observations apply to varia- 
tions of the heavy quark masses m c ,b, which, in principle, 
need to be considered as well, see Fig. [5] 



-0.001 



-0.002 



luj < 0.35, p Te > 1 GeV 
1.2<lT| K2.2,p > 1 GeV 




DSSV 




m e(1 [GeV] 



17: Impact of variations of the scales [if = fi r = 

+ [(Pt) 2 + (Pt) 2 ]/ 2 ) 1/2 on tne double-spin asymmetry 
A e j^ L for electron-muon correlations at NLO accuracy. The up- 
per panel shows A^ L computed with the DSSV parton den- 
sities for £ = 1 (solid line) and in the range 1/2 < £ < 2 
(shaded band). The lower panel gives the relative deviation 
of A^ L for £ = 1/2, 2 with respect to the result obtained for 
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Finally, we estimate the range of momentum frac- 
tions x at which the (un)polarized parton densities in 
Eq. are predominantly probed for the different single- 
inclusive and lepton-lepton correlation observables dis- 
cussed in this paper. Figure [T8l shows the corresponding 
cross sections differential in x. In each case, we have in- 
tegrated over the angular acceptance for detecting elec- 
trons and/or muons with the PHENIX experiment at 
RHIC as well as over all transverse momenta of the de- 
cay leptons. As before, an additional cut p^ > 1 GeV 
is imposed for lepton-lepton correlations. All x distri- 
butions are normalized to the respective integrated cross 
section (A)tTi nt . 

Since small transverse momenta probe also the small- 
est possible x values, Fig. [TS] gives a rough idea of the 
lowest possible momentum fractions accessible in heavy 
quark hadroproduction at RHIC with a c.m.s. energy of 
= 200 GeV. As can be seen, single-inclusive observ- 
ables, shown in the panels on the left, receive contribu- 
tions from x values as low as 10 -3 , but the majority of 
events has (x) ~ 0.01. Both, electron-muon and muon- 
muon correlations, displayed on the right-hand-side of 
Fig. 1181 probe on average larger values of x. In partic- 
ular, the latter observable can be used to study parton 
densities at (x) ~ 0.1 or higher. The entire suite of pos- 
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FIG. 18: Typical range of momentum fractions x at which 
the parton densities in Eq. ([7J are probed for the integrated 
unpolarized (upper row) and polarized (lower row) cross sec- 
tions for leptons from heavy flavor decays at RHIC shown 
in Figs. [2][5] The left panels give the results for the single- 
inclusive electron (solid lines) and muon (dashed lines) cross 
sections. Results for the electron-muon (dashed lines) and 
muon-muon (solid lines) correlations are presented in the right 
panels. 



identified through their semi-leptonic decays into elec- 
trons or muons. Such measurements can be done once 
sufficient statistics has been collected in polarized proton- 
proton collisions at RHIC and do not require the com- 
pletion of vertex detector upgrades. Decay lepton cor- 
relations turned out to be particularly suited probes for 
the spin structure of the nucleon, and experimentally rel- 
evant double-spin asymmetries of about one percent can 
be expected even for present-day gluon polarizations of 
rather moderate size. Contrary to naive expectations 
based on unpolarized results, gluon-gluon fusion is not 
necessarily the dominant channel for heavy flavor produc- 
tion in longitudinally proton-proton collisions at RHIC. 

Next-to-leading order QCD corrections are in gen- 
eral more pronounced for unpolarized than for polarized 
heavy flavor production cross sections, such that they do 
not cancel in the ratio defining double-spin asymmetries. 
Also, theoretical uncertainties estimated by varying fac- 
torization and renormalization scales are usually smaller 
for spin-dependent observables. 

The technical methods and their implementation into 
a parton-level Monte Carlo program can be straightfor- 
wardly extended to deal also with the spin-dependent 
photoproduction of heavy quarks. This will allow one 
to analyze already existing data for charm production 
consistently at NLO accuracy within future global QCD 
analyses of polarized parton densities. In addition, one 
can assess the physics impact of heavy flavor distributions 
and correlations obtained at a possible future polarized 
lepton- nucleon collider like the EIC. 



sible observables related to heavy flavor hadroproduction 
at RHIC can cover a wide range in x and has the poten- 
tial to provide novel information on the spin structure of 
the nucleon and the applicability of perturbative QCD in 
polarized hard scattering processes. 

IV. SUMMARY AND OUTLOOK 

We have presented a flexible parton-level Monte Carlo 
program to compute heavy flavor distributions and cor- 
relations at NLO accuracy in longitudinally polarized pp 
collisions. Experimental acceptance cuts, the hadroniza- 
tion of the produced heavy quark pair into D and B 
mesons, and their subsequent semi-leptonic decays can 
be included in phenomenological applications. 

Heavy flavor hadroproduction receives its importance 
for the field of spin physics from its partonic hard scat- 
tering processes, which differ from their counterparts for 
light hadron and jet production utilized in global QCD 
analyses so far. Upcoming measurements of heavy flavor 
production at RHIC will further our current understand- 
ing of factorization in the presence of spin and test the 
notion of universality for polarized parton densities. 

We have performed a comprehensive phenomenologi- 
cal study of various observables where heavy quarks are 
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Appendix 

In the Appendix we collect for completeness some ad- 
ditional details of the calculation and some explicit ex- 
pressions which were omitted in [171 ] but might be useful 
for the reader. 

First of all, we recall the LO partonic cross sections 
for open heavy flavor hadroproduction, which are needed 
for the factorization of collincar singularities. They also 
emerge in the soft gluon limit of the NLO 2 — > 3 real emis- 
sion contributions. The spin-dependent, color-averaged 
matrix elements squared for the tree-level processes in 
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([6]) read in d = 4 — 2e dimensions: 



A|Af 33 | 2 



A|Af 99 -| 2 



2Cp — Ca 



2t\u\ 



AB, 



QED 



\2 Cf 



where 



AB 



QED = 



AA 



QED 




(21) 
(22) 

(23) 
(24) 



Here, Nc denotes the number of colors, Ca = Nc, and 
Cf = (Nn — l)/(2iVcr). Contrary to the unpolarized case 
[la . Ha |, ABq ED receives no 0(e) contributions. The 
Mandelstam variables used in Eqs. (|2"Tj) - are defined 
by 

•5 = (Pi +P2) 2 , 

h = (Pi - ki f - m Q , 

ui = (pi - k 2 f - m 2 Q , (25) 

where s + 1\ + u\ = 0. p\£ are the momenta of the in- 
coming partons, and k\ and fc 2 are the momenta of the 
produced heavy quark and antiquark, respectively. To- 
gether with the appropriate two-body phase-space factor 
ifPS2 in d dimensions UM . 



mass parameter p is introduced to keep the strong cou- 
pling dimcnsionlcss in d dimensions. In the limit e — ► 0, 
F £ reduces to l/(167rs 2 ). In the c.m.s. frame of the in- 
coming partons, dPS2 and dAa a b can be conveniently ex- 
pressed in terms of the scattering angle between p\ and 
k\ by using the relation 



ti = --(l-/3cos0O, (29) 
where /3 2 = 1 — 4toq/s = 1 — p. The corresponding 
unpolarized Born cross sections in d dimensions can be 
found in Ref. [ll|. 

Next, we give explicit expressions for the p and u) pre- 
scriptions introduced in Eq. (|11[) to regularize soft and 
collinear regions of phase-space. The distributions are 
defined as follows 12311 : 



h(x) 



dx = 



1 h(x)-h(l) 



\1 x ) - J p 1 x 



dx, 



[h( X )-h(l)] l *f^±dx, 



X 

%) ( _L. ) 

i-w \ l -yJ w Ji-w l ~y 

-l+w , 1 X 

1 \i + 



fe(y)-M-i) 
1 + 1/ 



(30) 



dPS 2 



2tt 



[(47r) 2 -r(i- £ )]- 



tlUl — VTIqS 



x S(s + t\ + ui)dtidui , (26) 
the <?<? and qq Born cross sections can be written as 
d 2 Aa ab 



dt\dui 



= F s A\MabrS(s + t 1 +ui), 



where 



F^4[(47r) 2 -r(i- £ ) 



9 

jJLS 



(27) 



(28) 



collects all phase-space factors given in Eq. (|26[) and the 
flux factor l/(2s) included in Eqs. (EH) and The 



/i denotes an arbitrary "test function" , which is suffi- 
ciently regular in the limits x — > 1 and j/ — ► ±1. In 
a numerical implementation, the distributions (|30[) need 
to be used only if x and ?/ are sampled in the ranges 
p < x < 1 and 1 — w < y < 1 or — 1 < y < — 1 + u>, 
respectively, in the Monte Carlo integration. 

Following Ref. (HI, the soft, x — ► 1 limit of the 
2 — > 3 matrix elements squared is derived by applying the 
eikonal approximation, i.e., by taking the limit when the 
momentum of the emitted gluon gets soft. The obtained 
soft matrix elements squared for qq and gg scattering 
have the same form as in the unpolarized case [TEj but 
with the Born results replaced by their polarized counter- 
parts AAqed and AB QE d, given in Eq. ([24J) and 1(53]). 
respectively. The expressions can be straightforwardly 



integrated over x and y in (jHJ), and we obtain for dAa, 



0) 
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dAai 3 ) 



dtidui 



dt\dui 



Q 



2Li 2 1 



t-iU 



8Li 2 1 



I'Ul 



N c 



F £ G e a s (3~ AAqed \ -In 



2Li 2 1 - 



6Li 2 1 



t\Ul 

h 



tlUl 

^8Li 2 



8, h 
— in — 



2 sm Q 



^-hi 2 ^ 
e u\ t\U\ 



6Li 2 [ 1 - 

xti 



I In x In h In' 

Ml flUl 



g(2m| - s) 



[lnx-S^x)] j<5(s + + 

(31) 



6 In x In ■ 



In In x 



g(2m| - s) 

s(3 



S{h) U(a + *i+«i), (32) 



where we have used 

S(x) = ilnx+21nxln(l-x 2 )-ln 2 x 

+ Li 2 x 2 - C(2) , (33) 

G £ = Q4ire- e ^ E - ln ^ (1 - |c(2)e 2 ) f^f-) , 



and x = (1 — (3)/(l + (3). The dilogarithm function 
Li 2 (x) is defined as in Ref. [41| . and C(2) = 7r 2 /6 de- 
notes the Riemann Zeta function. G £ in (|34|) parame- 

I 



terizes the difference of the 2 — > 3 and 2^2 phase- 
space factors, the latter given by F e in Eq. (f2"5)l . As in 

Ref. [l5[ , we have split up the result for dAa^ into con- 
tributions from different color structures. The results for 
dAa^p and dAa^ A in Eq. d3TJ) and d32j) , respectively, 
agree with the corresponding unpolarized expressions in 
Ref. [TH after replacing AAqed by Aqed- With the 
help of (EH), dAa^ /dhdui can be easily transformed 
into dAa^q / d cos 9\ used in Eq. (| 13|) . 

Likewise, we obtain for the different color factors con- 
tributing to dAa\ S g : 



dAo 



gg,QED 



dt\dui 



dAa 



99,-4 



dt\dui 



dAa 



994 



dt\du\ 



4G 2 



N, 



2toq — s 
s/3 



r 2 



iln^ + V^ 

e iiUi 2 tiui 



In 2 — - In 2 x) + Li 2 ( 1 
ui \ 



[In x - S(x)] | <S(s + ii + ui) , 
_ 2fiui ' 

7 '<r 



In — [ -- +lnx) -Li 2 (l - — ) + Li 2 (l - 

ui V e / V uixy V fix 



2toq — s 
s/3 



2*itt 



- lnx — S(x) 



iV 2 



^(s + fi+ui), 
2 2tiui / 1 



2(2to^ 



s(3 



hui 



hi. 



1 + 



s 

hui 



In 



2 
Ul 



In 2 x 



S(x) 



>£(s + ti + ui) 



(35) 



(36) 



(37) 



Note that we have chosen a slightly different way of orga- 



nizing the above results according to their color structure 
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than in Ref. [Hj], but the sum of Eqs. ([35|) - ([37f agrees 
with their expression after replacing ABqed by its un- 
polarized counterpart Bqed [H[. The expressions for 
the virtual corrections dAa 1 ^ and the finite contribu- 
tions dAa^ in (fT"5|) are too long to be presented here. 
They are available upon request. 

Turning to the collincar, y — > ±1 limit of the 2 — > 
3 processes in ([3]), we give explicit expressions for the 

functions A/^^ appearing in Eq. (|14p . They read 



A/io +) (-, 



A/£ +) (Mi) 



h) = 327ra s s(l — .x 



AP- 



x A|M 99 | 2 | Pl 
327ra s s(l — .x 
x A|M 99 | 2 | P2 
327ra s s(l — x 
x A|Af 99 | 2 | Pl 

327TQ! s s(l — X 

x A|M W -| 2 | P2 

327TQ! s s(l — X 

x A|M, ? | 2 | P1 

327TQ! s s(l — X 

x A|M W -| 2 \ P2 ^ XP2 AP qq (x 



VXPI AP ff g(x 

>xp2 AP gg (x 
>a;pi AP ffg (a; 
■£cp 2 AP gg (x 



-.i'/j 



! AP M (x 



(38) 
(39) 
(40) 
(41) 
(42) 
(43) 



where the Born matrix elements squared in Eqs. (|21|) and 
(|2"2")l are to be evaluated with "shifted" kinematics. This 
is due to the collinear emission off one of the incoming 



partons such that only a fraction x of their original mo- 
menta pi t 2 is available in the subsequent hard scattering. 
The d = 4 — 2e dimensional LO polarized splitting func- 
tions APij in Eqs. (f3"B")) - f4"3"|) can be found in [30| and 
read: 



AP gg (x) 



AP gg (x) 
AP gg (x) 
AP 53 (x) 



1 - x + 3e(l - x) 



+ 5(1- x) ( ^£ + 21og/? 



I[2x-l-2e(l-x)], 

C F [2 - x + 2e(l - a;)] , 
1 

2x + 1 + 2e(l 



2CU 



'Ml -,■) ( ^+e^ + 21og/3j . (44) 



with /3q = llCyi/3 — 2n//3, fij as the number of ac- 
tive flavors, and where we have expressed the standard 
1/(1 — x) + distributions in AP qq and AP sg by the corre- 
sponding p-prescriptions defined in (|30[) . This amounts 
to introducing an additional log/3 term in the soft 6(1— x) 
parts of AP 99 and AP gs in ([13]). Note that in Eqs. 



(|43j) . contributions proportional to 6(1 — x) do not show 
up as they are already included as 1/e 2 poles in the soft 
cross sections listed in Eqs. (|3"Tj) - (l3"Tl) . In the factoriza- 



tion counter term dAa^ b in (|16p . only four-dimensional 
splitting functions are needed, i.e., s — > in (|44|) . 
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